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ABSTRACT 


Expected  information  gain  aa  a  result  of  life  testing  n  units  for  time 
t  is  calculated  for  the  time  transformed  exponential  model  and  a  utility 
function  based  on  entropy.  We  show  that  the  expected  information  gain  is 
concave  increasing  in  n  and  a  transform  of  the  test  time  t  .  A  computer 
program  for  calculating  expected  entropy  for  the  Weibull  distribution  model 
is  given.  TJ*is  may  provide  practical  guidance  in  designing  life  test  ex¬ 
periments  .  - - 
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EXPECTED  INFORMATION  FROM  A  LIFE  TEST  EXPERIMENT 


Richard  E.  Barlow  and  Jaw  Huan  Hsiung 


1.  INTRODUCTION 


In  considering  a  life  test  experiment,  two  questions  to  be  answered 


How  many  items  should  be  tested? 


How  long  should  we  be  prepared  to  wait  before  analysing  the  data? 


As  Llndley  (1956)  pointed  out,  "the  object  of  experimentation  is  (often) 
not  to  reach  decisions,  but  rather  to  gain  knowledge  about  the  world.” 
Hence,  we  do  not  consider  the  cost  of  experimentation  directly  in  a  con¬ 
ventional  decision  analysis  approach  to  the  solution  of  our  problems. 
Instead,  we  consider  the  Influence  of  sample  else  and  test  time  on  various 
measures  of  expected  information  to  be  gained. 

Since  the  objective  of  testing  is  to  gain  information  about  life 
times  of  similar  items,  we  need  to  determine  how  our  expected  measure  of 
Information  to  be  gained  depends  on  sample  sice  as  well  as  test  time. 


By  information ,  we  mean  anything  which  changes  our  probability  distribu¬ 
tion  about  unknown  quantities.  To  measure  this  change  we  use  a  utility 
function,  u(A,d(D))  ,  where  A  is  the  unknown  life  distribution  para¬ 
meter  of  interest  and  the  decision  taken,  d(D)  ,  based  on  observed  data 
D  ,  will  (in  this  paper)  usually  be  identified  with  the  posterior  mean 
or  the  posterior  density.  The  expected  gain  in  information  based  on  n 
observations  can  than  ba  meaaured  by 
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g(n)  -  E  max  J*  u(X,d)ir(X  |  D,n)dxj  -  max  J  u(X,d)*(X)dX  (1.0) 


where  A  is  the  parameter  space,  n  is  a  prior  density  for  X  and 
ir(X  |  D)  is  the  posterior  density  for  X  given  data  D  and  d  belongs 
to  some  appropriate  decision  space.  Raiffa  and  Schlaifer  (1961)  call 
(1.0)  the  expected  value  of  sample  Information.  The  expression  is  easily 
seen  to  be  nonnegative.  This  idea  of  measuring  expected  information  aa 
expected  utility  has  been  discussed  by  DeGroot  [(1970),  pp.  429-433]  and 
more  recently  by  Bernardo  (1979). 

To  illustrate  ideas,  first  consider  a  non- life  test  situation  where 

n  normally  distributed  measurements  are  to  be  made.  Suppose  our  uncer- 

2  2 
tainty  about  measurement  X  given  0  and  o  ia  measured  by  a  H(0,o  ) 

2 

distribution.  For  convenience,  suppose  o  is  known  but  6  is  unknown 
so  that  we  wish  to  learn  about  6  .  Let  our  prior  uncertainty  for  6  be 
measured  by  a  distribution.  Let  x^^*  ••.,  x  be  n  in¬ 

dependent  (given  6)  observations  so  that 


x  -  (Xj  +  ...  +  x^/n 


given  8  has  a  distribution  while  9  given  x  has  a 

fl{y(x),Tg)  distribution  with  mean 


y(x)  ■  (1  -  w)8o  +  wx 


and  variance 


ft 


where  v  * 


2 


.  To  measure  information  gained  as  a  result  of  these 


n 


Y2  + 


measurements. 


a 

n 

let  our  utility  function  be 


u(8,d) 


1 

0 


if  |6  d}  <  e 
otherwise  . 


Then 


(1.1) 


m 
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u(e,d)*(0  |  x,n)d8 


-  Maximal  F(|«  -  d|  <  t  |  x,n] 
d 

•  P[|«  -  u(x)|  <  c  |  x,n] 


so  that  d  •  p(x)  Is  our  optimum  "decision"  in  this  case.  Therefore, 
if  we  take  n  measurements,  our  expected  utility  will  be 


I 


u(6,ii(x))v(6  |  x,n)de| 


*{p(|e  -  u(i)|  <  e  I  i,n]} 


(1.2) 


This  is  our  expected  posterior  probability  that,  after  n  measurements 
are  made,  6  #111  be  within  e  of  the  posterior  mean  y(x)  .  ♦  is  the 

cumulative  M(0,1)  distribution.  One  way  to  determine  n  irrespective 
of  cost  considerations  is  to  specify  a  probability  p  and  raquire  that 
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then 
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and 


where  [• ]+  denotes  greater  nonnegative  Integer  In  the  quantity  within 
brackets.  If  we  let  y  •  •  >  then  we  have  the  non-Bccyeeian  solution 


Pros  (1.1)  we  see  that  our  expected  measured  gain  in  information 
baaed  on  n  observations  will  be 


gin) 


-4/ 


E{  J  u(6,ti(x))v(0  |  x,n)d0 
® 


u(0,0  )«(0)do  (1.3) 

o 


where  0Q  is  the  prior  mean  and  y (x)  is  the  posterior  mean.  It  is 
easy  to  verify  from  (1.2)  that  g(n)  is  concave  increasing  in  n  so 
that  marginal  gain  is  decreasing  in  sample  size  n  . 

There  are  several  reasons  why  s  utility  function  such  as  (1.1)  and 
the  expected  measured  gain  (1.3)  based  on  (1.1)  might  not  serve  as  an 
adequate  measure  of  information  gained. 


[log  ff(9)}*(e)iie  a 


firat  introduced  end  ataadiad  by  Lindlay  (1956) .  It  is  tha  negative 
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* 


i 
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where  p(D,6)  is  the  joint  density  of  data  D  and  parameter  6  while 
p(D)  *  J  p(D,0)d6  .  If  0  *  g(w)  is  a  1-1  map  of  (5)  onto  (5)  and  J 

is  the  Jacobian  of  the  transformation,  then 

*W  "  JJp<b>«<w»IjI  lo«  [pwiu(l))]dD<)v 

.//p(D.g(«))|j|  log 

■//  P*(D-“)  lo*  [-ptn)^)]^ 

where  p*  is  the  joint  density  of  D  and  w  .  Hence  (1.4)  is  invariant 
under  1-1  transformations  of  the  parameter  space. 

We  will  show  for  the  time- transformed  exponential  life  distribution 
model  and  the  utility  function 

u(0,p(*))  -  log  p(0) 


where 


# 


WHY 

p(0 


[log  p(6)]v(0  |  D)d0 


[log  ir<0  |  D)]ir<0  |  D)d0 


that  our  expected  measured  gain  in  information  is  concave  increasing  in 
both  sample  size  n  and  a  transform  of  the  test  time  t  .  Methods  for 
calculating  expected  information  with  respect  to  a  Weibull  life  distribu¬ 
tion  model  are  discussed. 


7 


INFORMATION  FROM  A  LIFE  TEST  EXPERIMENT 
-Ro(x) 

Let  Fq(x)  *  e  be  a  specified  absolutely  continuous  life  dis¬ 

tribution  with  hazard  function  R  and  failure  rate  r  (x)  *  ~r~  R  (x)  . 

o  o  ox  o 

Consider  the  life  distribution  model 


-XR Ax) 

F(x  {  X)  -  e  °  (2.1! 

where  X  is  the  unknown  "proportional  hazard"  but  RQ  is  specified. 
(2.1)  is  called  the  time-transformed  exponential  life  distribution  model 
Suppose  n  similar  units  are  put  on  life  test  for  the  time  interval 
[0,t]  and  we  judge  the  model  (2.1)  to  be  an  appropriate  description  of 
our  uncertainty  concerning  the  life  length.  If  we  observe  k  failures 
with  lifetimes  x^fX2»  ...»  x^  and  n  -  k  survivors  in  [0,t]  ,  then 
the  likelihood  is 


L(X  I  *i»X2*  *k»t)  “  (k)^^  W] 

[’x[ii  R°<Xi> + <n‘  k)Ro(t>]]  * 


(2.2: 


exp 


k 

Clearly  k  and  s  »  T  R  (x  )  +  (n  -  k)R  (t)  together  constitute 

i-1  0  1  0 

a  sufficient  statistic  for  X  .  For  some  results,  we  will  use  the  prior 
density 


*(X) 


by*.-** 

ru) 


X,a,b  >  0 


(2.3! 


and  posterior  density 


8 


ir(X  |  k,s) 


(b  +  s)*^*4*"1  -(b+s)X 

rcT+T)  c 


2.1  A  Measure  of  Information  Based  on  Entropy 

Lindley  (1956)  introduced  the  following  measure  of  expected  informa¬ 
tion  gain  as  a  result  of  performing  an  experiment  E  resulting  in  data  D 


I(E,ir(X))  -  E  J  [log  ir(X  |  D)]ir(X  |  D)dX  -  J  [log  ir(X) Jir(X)dX  (2.4) 
A  A 

where  the  expectation  operator,  E  ,  is  with  respect  to  the  unconditional 
distribution  of  the  data  D  .  Bernardo  (1979)  pointed  out  the  connection 
with  expected  utility  where  the  utility  function 

u(X,ir(X  |  D))  ■  log  ir(X  |  D)  (2.5) 

depends  on  X  and  the  decision  variable  is  ir(X  |  D)  ,  the  posterior 

density  at  X  .  The  entropy  is  -/  [log  ir(X  |  D)Jir(X  |  D)dX  and  (2.4) 

A 

is  the  negative  expected  change  in  entropy  as  a  result  of  performing  E  . 

Information  measures  based  on  (2.4)  are  dimensionless  and  as  such 
may  be  difficult  to  interpret.  However,  (2.4)  does  provide  a  vay  of 
ordering  proposed  experiments  by  assigning  information  values  which  are 
invariant  under  1-1  transformations  of  the  parameter  space.  For  example, 
suppose  we  life  test  n  units  for  time  t  and  use  the  Weibull  life 
distribution  model 


P(X  >  x  |  a,X) 


(2.6) 


where  a  is  known  but  X  unknown.  Also  let  the  prior  for  X  be 


ir(X  |  A,B)  -  BAXA"1e"BX/r(A)  . 


Figure  2.1  and  2.2  are  example  graphs  of  expected  information  versus 
sample  size  and  test  time  respectively.  For  example,  from  the  graphs, 
we  can  see  that  testing  3  units  for  3  years  results  in  the  same  infor¬ 
mation  as  testing  10  units  for  about  0.75  years.  Thus, we  have  a  means 
of  comparing  experiments.  Information  values  can  be  related  to  specified 
experiments. 

The  parameters  of  the  gamma  prior  used  (A  and  B)  were  originally 
specified  based  on  the  pressure  vessel  data  analyzed  in  Barlow,  To land 
and  Freeman  (1979).  The  shape  parameter  a  -  1.5  was  used.  The  graphs 
show  that  for  these  parameter  values  (A  and  B)  there  Is  little  to  be 
gained  by  testing  more  than  3  years. 

In  order  to  obtain  our  main  results  we  define  the  experiment  E 
as  a  quadruple  {0,8, A, P}  ,  where  P  is  the  space  of  observations  x 
of  the  random  vector  X  ,  8  is  the  a-field  of  the  subsets  of  P  ,  the 
probability  measure  (or  density  of  X  belongs  to  a  family  P  indexed 
by  a  parameter  X  e  A  .  Suppose  that  the  observation  x  in  our  ex¬ 
periment  E  consists  of  a  pair  of  observations  x^fX2  *  *****  **' 
p  ■  P^  x  P2  .  Let  8^  be  the  o-field  over  P^  induced  from  8  by 
the  transformation  x^  -  x^(X)  and  let  P^  be  the  set  of  probability 
measures  on  81  (i  ■  1,2)  .  Then  E^  -  {P^,8^,A,P^J  (i  -  1,2)  are 
two  experiments.  Denote  the  sum  of  the  experiments  E^  and  E 2  » 
by  E  ■  (E^E^  .  Now  we  consider  a  related  experiment  E^x^)  ■ 
{D2,B2,A,P2(x1))  ,  where  P2(x2)  is  the  set  of  probability  measures 


of  Xj  conditional  on  x^  .  Mow  consider  the  expected  information  for 
E2  were  we  to  know  the  observation  x^  from  performing  : 


I  xl^  "  V  j  l1®*  *(*  I  *2»Jtl^1r^A  I  Xj.XjMl 
*  A 

(2.7) 

“  J  (log  r(A  |  XjJJirCX  |  Xj)dX  . 

A 

Since  »(X  |  Xj^)  ia  the  posterior  density  of  X  after  Xj  has  been 
observed*  I(E2<x1)  »v(X  |  Xj))  is  the  measure  of  expected  Informal ion 
gain  to  be  provided  by  our  Observation  *2  after  E^  has  been  per¬ 
formed  and  xx  observed.  I(E2  |  E^)  -  E^[I(E2(x1),»(X  |  x^)]  , 

the  average  of  I(Ej (*j>.»(*  }  Xj))  over  Xj,  ,  la  defined  to  be  the 
average  information  to  be  provided  by  Ej  after  ^  has  been  performed. 
From  now  on  we  shall  often  denote  the  expected  information  by  1(E) 
whan  the  particular  prior  distribution  deee  not  have  to  be  stressed. 

This  measure  of  information  has  the  following  properties: 

1.  1(E)  >  0  . 

2.  1(E2  |  Ex)  >  0  . 

3.  1(E1)  +  I(E2  |  Ej)  -  1(E)  whera  E  -  (EpEj)  . 

4.  If  xx  is  sufficient  for  X  ,  then  I(E2)  -  1(E)  . 


_  ... 


;.-.’.-‘.,.-.v  .v'.v-v.v  v 
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5.  If  p(xiyx2  |  X)  ■  p(x2  |  X)p(x2  |  X)  ,  i.e.,  x2  and  x2  are  In¬ 
dependent  When  X  Is  given,  then  I(E2  |  E^)  <_  I(E2>  . 

6.  Let  ■  E^  be  any  experlaent  and  let  E^E^  ...  be  independent 

identical  experinsnts.  Let  E^  •  (E^,E2>  and  generally  E^  ■ 
(En’E(n-i))  *  Th*n  *  concave  increaaing  function  of  n  . 

See  Lindley  (1956)  for  proof a  of  the  above  properties. 

Let  E_  .  „  be  the  experlaent  wherein  n  units  aged  t.  and  with 

identical  life  distributions  arc  put  on  life  test  to  age  t2  (t2  >  t,)  . 
Assuae  statistical  independence  aaong  the  n  units  conditional  on  X  . 

Thaoran  2.2 

_  “XR  (x) 

For  the  tlaa-transforaad  exponential  nodal,  F(x  |  X)  •  a  , 

I<E«,0.t>  *•««»«  !«««*«•  *•  »  «d  =»««.  l«n»la(  n. 

RQ(t)  .  The  prior  density  *(X)  is  arbitrary. 

Proof: 

Since  1  (•)  is  known  and  continuous,  T  ■  R  (X)  la  exponentially 
O  P. 

distributed  with  paraaater  X  .  Therefore,  perfomlng  an  experlaent  for 
a  period  [0,t]  under  the  tlae-traaaforaed  exponential  nodal  ia  the  sane 
aa  parforaJLng  an  experlaent  for  a  tiaa  period  [0,RQ(t))  under  the  ex¬ 
ponential  aodel.  Let  *(X)  be  the  prior  in  both  caaea.  Hence, the  ana e- 
urea  of  expected  laforaatlon  gain  provided  by  these  two  cxperlaants  are 
the  sans.  It  la  therefore  euffldant  to  prove  thla  theorea  for  the  ex¬ 
ponential  nodal. 
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How  define 


:<0)  . 


n,0,t-A  * 


E«>  -  E 


n-k0,t-A,t  *iven  ko  failure8  in  I0,t-Aj  , 


:<2) 


a-kg-k^.t.t+A  given  Kq  failures  in  [0,t-A] 
and  k^  failures  in  { t-A , t ]  . 


+ 


r(0) 


+ 


:(1) 


+ 


:(2) 


t-A 


t4A 


n,0,t-A 


En-k0,t-A,t  En-k0-kltt,t4A 


Than 


I(En.0.t/<A»  “  I(En,0,t-A»1r(A)) 

I(E<»  |  E<°>)  -  |  D0» 


using  tha  memoryless  proparty  of  tbs  exponential  and  DQ  ■  (kg  ,  total 
tine  on  teat  In  (0,t-A))  ,  the  sufficient  statistic  for  X  .  Similarly 


X(E<»  |  £<».£<»)  -  I  VV> 


where  •  (k^  ,  total  tlae  on  test  in  (t-A,t)) .  To  show  concavity  and 
the  Increasing  property  we  need  only  show 


0  <  1(B(2)  |  E(1),E(0))  <  I(E(l)  |  E(0))  . 


-v-v- .  •  >v  v  .'-v-  .-.  .v-v-v  v-v 

■■  - 8 atfadi -----  i -- 


v  r;  t  ,  r.  v; 
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By  definition  and  property  (5) ,  we  have 


1  do*bi»> 

■  I<En-k0,t,t+4,,,<1  I  V  I  ED-k0,t-4,t,1'(>  I  V* 
i  1<£.-lc0,t,t«-’<‘  I  D0»  • 


By  using  the  aenoryless  property  of  the  exponential  end  using  the  prior 
ir(A  |  Dq)  on  the  parameter  space,  we  have  KE^^  (,t,t4A,ir^X  I  V>  “ 

I(En-k0.t-A,t,1,(X  I  D0))  •  Th#r*fore» 

ED1(I(En-k0,t,t+A,w(X  I  D0*V))-I(En-k0,t-A,t#ir(X  I  V*  * 

But  information  Is  Increasing  In  sanple  sice,  so  that  if  fc^  is  the 
(random)  nuaber  of  failures  in  E  then 

I  D0-V))iI<E»-k0,«-4.f*<X  I  D0»  ' 


Mow  take  the  expectation  with  respect  to  Dn  .  Then 


0<X(E«>  |  E(W.E<0>)  <  i(£<»  |  E<°>>  . 

Hence,  I(E_  A  „)  Is  concave  Increasing  In  t  which  coapletes  the  proof. 
n,u,t 

2.2  A  Cowouterlzed  Method  for  Calculating  Entropy  in  the  Case  of  a 


P (Ilf seise  >  z)  •  exp  (-Ax)  ,  where  x  >  0,a  >  0,A  >  0  . 

This  is  Che  Welbull  distribution  survive!  probability.  Assume  a  is 
known.  From  (2.2),  the  likelihood  function  of  A  is 

L(x  I  *1 . V*>  *  (k)lk“k[1I1  *l]  “p  Wjj  *i  *  <n  '  k)t°]|  • 

k 

The  peir  k  end  s  »  T  x,  ♦  (n  -  k)t  constitute  e  sufficient  stetlstic 

i-1  x 

for  A  .  Let  K  end  S  be  the  rendon  quentitlee  corresponding  to  the 
number  of  fellures  end  totel  time  on  test,  respectively.  Bartholomew 
(1963)  hes  obtained  the  joint  density  of  K  end  S  given  A  es  follows: 

?(K-k,S-.  |  X)  -  (1)  .*X*  !  [0..- t“<i.-k+i)]jk  1 


5  Dk(s)Ake~A8 


where 

V>  -  (9  7k^T)T  X  (i)w»l[—  1°*»  -  *“<»  -  k  +  1)1  r*  • 

The  probebllity  of  observing  no  feilure  in  [0,t]  is 

P[K  -  0  ,  S  -  nt®  |  A]  -  e'Ant°  s  p(0,nt“  |  A) 
where  's'  means  definition. 

Assume  (2.3)  es  the  prior  density  of  A  .  Using  Equation  (10)  in 


Lindley  (1956),  we  have 


v*a>. 


nt°  » 


where 


-  i  if  d,ax 


•■0  A«0 


nt  • 

■iff  p(k,o  |  A)r(A)  log  p(k,8  |  A)dsdX 

kaO  J-  V- 


•■0  x«o 


nt 

n  r 

-If  P<k. 

k-0  J. 


•)  log  p(k,s)d8 


m 

-  J  p(0,nt“  |  X)ir(X)  log  p(0,nt“  |  X)dX 


nt°  - 


iff  P(k»*  |  X)r(X)  log  p(k,«  |  X)dsdX 

k-1  A*  * 


8-0  X-0 


n  nt 

-  p(0,nt°)  log  p(0,nta)  -if  p(k,»)  log  p(k,s)ds 

**  8—0 

nt° 

S  AQ  -  Bq  ♦  J  [Aj^(s)  -  Bk(8)]di  , 


m 

A0  -  J  p(0,nt“  |  X)w(X)  log  p(0,nta  |  X)dX 


m 


XOt°  (“Xnt°)dX 


(b  ♦  nta) 


•+1  * 
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»Q  -  p(0,nt°)  log  p(0,nt°) 


Also 


I  ?  -Xnt°  b*X-1.-bX  .1  .  If  -Xnta  b4X*-1«-|,X  .. 
"  J  4  ~  ru~  dx  lo»  J  *  r(.)  dx 

|X«0  J  |x«o 

(b  +  nta)  |_ (b  +  nta)  j 

m 

v>*  J  p(k,s  |  X)ir(X)  log  p(k»s  |  X)dX  , 


Bk(»)  “  P(k,»)  log  p(k,s)  . 


/V  la 

Dk(s)XKe  AS  b-  [log  Dk(s)  +  k  log  X  -  Xs]dX 


a. a+k-1 - (b+») X 


J.  ■ 

- ru 


«.«+k-l  - (b+s)X 


/.  a.*+k-l  -i 

4 - r5T 


log  XdX 


/h«i*+k  -(b+s)X 

- dx 

0 


r(«  +  k)b*Dfc(«)  log  D.  (•) 

- k - __k - +  A.’ 


ruHb  ♦  •)' 


•Dk(«)b*T(»  +  k  +  1) 


ruxb  ♦  •) 


a+k+1 


t-  ■  y  ■ 


|f»|&  ^nj 

♦  r j  .'  ’  ■  •> ' 

wSP^jHH 

|ggig 

onow  u>  o  o  n  o  nnn  h  nnnnnnnnn 


APPENDIX 


PROGRAM  INFO (INPUT .OUTPUT .TAPE5-INPUT .TAPE 6 -OUTPUT) 

THIS  PROGRAM  CALCULATES  THE  MEASURE  OF  INFORMATION  (EXPECTED 
ENTROPY)  OF  A  LIFE  TEST  EXPERIMENT  WHEREIN  THE  LIFE  DISTRIBUTION 
OF  THE  TESTING  UNIT  IS  WEIBULL  DISTRIBUTION  AND  PRIOR  IS  A  GAMMA 
DISTRIBUTION.  WE  USED  SUBROUTINE  ’GAUSS Q’  TO  EVALUATE  THE  INTEGRAL. 
’GAUSSQ'  APPLIES  GAUSSIAN  QUADRATUE  TECHNIQUES  TO  DO  THE  EVALUA¬ 
TION  OF  THE  INTEGRAL. 

THIS  PROGRAM  IS  GOOD  FOR  SAMPLE  SIZE  UP  TO  50. 

DIMENSION  BB(500).X(500). C(500).ENDPTS(2) 

COMMON  FACTA (55) .TALPHA.A.B.N.H 
READ  100.  N,  A.  ALPHA.  T.  B 
IF  (  N.EQ.O)  STOP 
IF  (N.  GT.  1)  GO  TO  3 

INFORMATION  CALCULATION  FOR  N-l 

TA-T**ALPHA 

V-B+TA 

VB-B/V 

VA— VB**A 

UA—  (A*TA*VA)/V 

UB-VA*ALOG(VA) 

Al-A+1 

UC- (1— VA) * (PS I ( Al) -ALOG (A) ) 

UD-A*VA* ( 1-VB+ALOG (VB) ) 

XINFO-UA-UB+UC+UD 
PRINT  150  ,N,T,  ALPHA.  A.B 
PRINT  300.XINF0 
GO  TO  1 

INFORMATION  CALCULATION  FOR  N  GREATER  THAN  1 
GENERATE  FACTORIAL  FROM  0  TO  N 

FACTA(1)-1. 

FACTA(2)-1. 

MN-N+1 
DO  5  1-3. KN 
K-I-l 

FACTA(I)-K*FACTA(K) 

CALCULATE  AO  AMD  BO 

TALPHA-T**ALPHA 
H-N^T ALPHA 
Q-B+H 

QQ1-  (B/Q)**A 
QQ2-(A*H) /Q 
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AD—  (QQ1)*QQ2 
B0-(QQ1)*AL0G(QQ1) 

PRINT  200 ,N,T, ALPHA, A,B 
CO- AO- BO 
PRINT  250,  CO 
KIND-1 
KPTS-0 
DO  25  1-2,3 
MM-90M 
C 

C  WHEN  N*TALPHA  IS  LARGE,  IT  IS  BETTER  TO  CHARGE  THE  VALUE  OF  MM, 

C  FOR  EXAMPLE  MM-90*I,  MM- 120* I . BUT  THE  VALUE  OF  MM 

C  CANNOT  GREATER  THAN  500. 

C 

CALL  GAUSSQ(KIND,MM,BALPHA,BETA,KPTS,ENDPTS,BB,X,C) 

C 

C  'GAUSSQ'  RETURNS  THE  NODES  X(I)  AND  WEIGHTS  C(I),  THEN  APPROXIMATES 

C  THE  INTEGRAL  BY  SUM  OF  C(I)*F(X(I))  (I  FROM  1  TO  N) . 

C 

D-0.0 

DO  20  J-l.MM 
20  D-D+C(J)*F(X(J)) 

PRINT  150,  D 
XINFO-CO+D 
PRINT  300.XINF0 
25  CONTINUE 

100  FORMAT(I3,F6.3,F5.2,F4.1,F5.2) 

150  FORMAT (//10X,9HINTEGRAL-,F21.14) 

200  FORMAT(///10X,2HN-,I3,10X,2HT-,F4.1,10X,6HALPHA-,F5.2,10X,2HA-,F6 

•3,10X,2HB-,F5.2) 

250  FORMAT (//10X,3HCO-,F21. 14) 

300  FORMAT ( //10X, 12H INFORMAT ION- ,F21 . 14) 

GO  TO  1 
END 


FUNCTION  F(S) 

C 

C  TO  USE  'GAUSSQ'  AN  INTEGRAL  (FROM  A  TO  B)  OF  F(X)  MUST  BE  BROUGHT 

C  TO  THE  STANDARD  INTEGRAL  FORM.  THIS  IS  DONE  BY  A  SUITABLE  CHANGE 

C  OF  VARIABLES,  FOR  EXAMPLE,  INTEGRAL  (FROM  A  TO  B)  OF  F(X)  EQUALS 

C  TO  (B-A) /2  TIMES  THE  INTEGRAL  (FROM  -1  TO  1)  OF  F(Z),  WHERE 

C  Z»(X+1) (B-A) /2+A. 

C 

DIMENSION  SUM(55) ,D(55) ,SUMM(55) 

COMMON  FACTA(55) ,TALPHA,A,B,N,H 
S-((S+l)*H/2. 

TIND-TALPHA* (N-l) 

IF  (S-TIND)  30,30,40 
30  D(l)-0. 

GO  TO  45 
40  D(l)-N 


DO  50  K-2.N 
SUM(K)-0 
DO  70  K-2,N 
DO  60  J«1,K 
TEST-TALPHA* (N-K+J-l) 

IF  (S.LE.TEST)  GO  TO  65 
TT1"(S-TEST) ** (K-l) 

TT2*FACTA(K+1) / (FACTA (J)*FACTA (K-J+2) ) 
TT3-(-l)**(J-l) 

SUM(K)  "SUH(K)  4TT1*TT2*TT3 

TT4-FACTA (N+l ) / (FACTA (K+l ) *FACTA (N-K+l) *FACTA (K) ) 
D(K)-TT4*SUM(K) 

DO  80  K-l.N 

TT5-(B/(S+B))**A 

TT  6-GAMMA  (A-HC)  /GAMMA  (A) 

TT7-D(K)/((S+B)**K) 

TT8-TT7*TT5*TT6 

X-A4K 

Y-PSI(X) 

TT9«K*Y- (S*X) / (S4B) 

TT10»TT9- ALOG (TT 5} -ALOG (TT6 ) 

SUMM(K) «TT8*TT10 
XSUMM-0. 

DO  90  K-1,N 
XSUMM-XSUMMf  SUMM  (K) 

F-(H*XSUMM)/2. 

RETURN 
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